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ABSTRACT 



Results on the statistical analyses of series of events published 
subsequent to the monograph by Cox and Lewis on this subject are surveyed 
Special emphasis is given to tests for renewal processes, a considerable 
amount being now known about the distributions of some of the test sta- 
tistics involved, and to testing the functional form of a trend in a 
nonhomogeneous Poisson process, as well as the point process model itself 
A survey of work in special processes such as cluster processes and 
doubly stochastic Poisson processes is also given. 



Prepared by: 



1. INTRODUCTION 



The object of this paper is to survey recent results in the statisti- 
cal analysis of univariate point processes (series of events). The survey 
is personal reflecting my own present interests, and is not comprehensive 
For convenience, I have taken '*recent** to mean anything published since 
Cox and Lewis (1966). I have also mentioned some topics and results 
which were omitted, for reasons of emphasis or ignorance, from that 

monograph. Finally, I point out areas where further work is required. 

} 

A survey of work in this field is not so difficult as a survey of, 
say, the theory of point processes because the advances in the statisti- 
cal analysis of point processes are, by comparison with the theory, few. 
This may reflect the relative difficulties of these areas, but that thought 
may be a personal bias. 

Some of the shortcomings of the Cox and Lewis monograph were, 
in hindsight, the following: 

(i) Not enough consideration of grouped data. In some cases where 

point events occur, such as epidemiology, recording limitations or the 
volume of data force one to work with numbers of events in fixed inter- 
vals whose width may or may not be controllable in advance. I have 
touched on this problem recently (Lewis, 1970)| separate spectral 
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analyses of the intervals and of the counting process (Bartlett, 1963) 
are, for example, not possible. One has to use a spectral analysis of 
the grouped counts (i. e. , number of events in a fixed time interval), 
which is very nonstandard in its distribution theory, but in which well 
known problems of aliasing arise. Results of Cox (1970) may be useful 
here, and Cox (1955) has considered some other aspects of the analysis 
of grouped events. 

Interesting examples of this type of problem are found in physics 
and optics. For example (Helstrom, 1964; Karp and Clark, 197 0) 
photon or other particle emissions are known from physical consider- 
atr>ns to be generated by a doubly stochastic Poisson process and it is 
required to determine parameter values of the driving process from 
counts of the number of photons emitted in successive periods. Here 
it is prohibitively costly to record exact times of occurrence. However, 
the recording interval can be determined in advance by the experimenter. 
Problems of grouped Poisson counts (McNeill, 1971) are also common. 

(ii) Very little emphasis was given to sequential methods. While this 
was to deliberate to save space, it is also true that most data I have come 
across is presented for analysis in toto. This may change as better re- 
cording methods are introduced and, hopefully, as statisticians are called 
in before the fact. There is still a problem in that simple sequential 
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methods are known only for Poisson processes (homogeneous or non- 
homogeneous) and also that rather unsmooth inhomogeneities occur in 
practice which make the application of formal sequential methods based 
on very definite assumptions quite hazardous. 

Less formal sequential methods are useful; in particular, an 
analysis of the data in successive sections is very useful, both to cut 
down on computation time in, say, spectral analyses (Lewis, 197 0) 
and to examine the time evolution of the process. 

As an example, consider a series of arrivals at an intensive 
care unit in a hospital. This data will be used for illustration through- 
out the paper; it was supplied by Dr. Barr, of the Oxford Region- 
al Hospital Board, England. The first section, consisting of 
n= 251 arrivals in t^ = 409 days (4 February, 1963 - 18 March, 1964), 
was analyzed in Chapter 1 of Cox and Lewis (196 6). Later on, the 
arrivals up to 6 February, 1968 were received. Three subsequent 
sections of length t^ = 409 days were taken from these later arrivals 
for comparison and their statistics, as well as those for the total 
record, are shown in Table 1. 

* These times -of -arrival were exceptionally well recorded. Of the 
1468 arrivals in the 1420-day period from 4 February, 1963 to 
6 February, 1968, only one time-of-day was not recorded. Nine tied 
arrivals occurred. Generally, recording times seemed to be at the 
five-minute intervals of the hour, although other times do occur. The 
data to 18 March, 1964 is given in Cox and Lewis (1966, p. 255); the 
rest in Lewis (1971). 
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Table 1. Arrivals at an intensive care unit (no ties) 



Period 1 

4 Feb '63- 
18 Mar '64 



Period 2 

19 Mar '64- 
2 May '65 



Period 3 

3 May '65- 
15 Jvme '66 



Period 4 

16 June '66- 
29 July '67 



Total 

4 Feb '63- 
6 Feb '68 



n 251 350 

*^'0 409 409 

U(days) 218.47 218.53 

U 1. 874 2.222 

maN/t^ 0.614 0.856 



372 316 1458 

409 409 1829 

183.69 197.20 954.25 

-3.399 -1.099 2. 874 

0.910 0.773 0.797 



n 

The statistic U = S t, /n, where the t. are the times to events and 

i*l 1 1 

n the number of events in the period of observation of length t^, is the 
centroid of the t.'s and is used to test (Cox and Lewis, 1966, p. 47) 
for the trend parameter p in a nonhomogeneous Poisson process with 



rate 



( 



X. (t) » exp(or + p t)* X exp(p t). 



( 1 . 1 ) 



In testing for p a 0 against p 0, <r is a nuisance parameter with 
sufficient statistic n and the test is based on the distribution of U, 
given n. The normalized statistic U in line three of Table 1 should 
be distributed approximately as a N(0, 1) variable. There was fairly 
strong evidence of a monotone trend at the end of the first observation 
period ( 7 percent level), but later on in the series, as in Section 3, 

there is a definite decreasing trend. However, the total arrival process 
gives fairly strong evidence of an increase in X (t), the value U » 2. 874 
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being significant at about a 2 percent level. For Period 1 and Period 2 combined, 

U was found to be 4.492 which is significant at a level much smaller than 1 percent. 
An actual sequential test of p = 0 against p 0 rejects the null hypo- 
thesis at a 1 percent level after 550 days. The sectioned analysis of 
four periods suggests a long cycle or quadratic term in the exponential 
trend (1. 1). The nonhomogeneity is confirmed very strongly by a dis- 
persion test (Cox and Lewis, 1966, p. 232) applied to the numbers of 
events in the four sections. This has value 26.74; its distribution is 
that of a chi-square variable with 3 degrees of freedom which has a 
0. 9999 point equal to 21. 11. 

A plot of the times -to-events t. against serial number i is 
given in Figure 1, and a plot of the average values of successive sets 
of twenty intervals between arrivals is shown in Figure 2. Again the 
long cycle or quadratic trend is quite evident graphically. We return 
to this example later in the paper. 

(iii) Sectioning brings up in some ways the analysis of replicated 
point processes, which again was not considered in detail by Cox 
and Lewis (1966). This replication can occur quite commonly in exper- 
imental situations, for instance in neurophysiology where experiments 
can be repeated many times. Here the signals are trains of very narrow 
spikes of apparently fixed height, so it is appropriate to analyze the 
times of occurrence of the spikes as a point process. Observation over 
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TIME TO ARRIVAL IN DAYS 
(EVERY 10) 

Figure 1, Arrival of patients at an intensive care unit. 

Complete record, 4 February 1963 to 8 February 1968, Number of arrivals (n) 1458 in 
period of t^ * 1829 days. Arrival number vs, time to arrival, m = n/t^ = 0. 7972. 



AVERAGE OF SUCCESSIVE GROUPS OF 20 
INTER-ARRIVAL TIMES 
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SERIAL NUMBER OF MIDDLE ARRIVAL IN GROUPS OF 20 INTER-ARRIVAL 
TIMES 



Figure 2. Arrival of patients at an intensive care unit. 

Averages of successive groups of 20 inter-arrival times. Nvimber of arrivals (n) 1458 
in period of t •= 1829 days. Overall average is 1. 264 days. Homogeneity of variance 
statistics » 318. 000 mean 71, a- = 11.92). 
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too long a period can be deceiving, because of physiological deterioration, 
and replication is sometimes preferable. Of course, care must still be 
taken that physiological or experimental conditions have not changed. 

Comparison of rates and trends, mainly in Poisson processes, 
was discussed in Cox and Lewis (1966). See also Qureishi (1964) for 
the comparison of rates in twoWeibull processes. General problems 
of multiple point processes (multivariate point processes) are discussed 
in Cox and Lewis (197 2) and Perkel, Ger stein and Moore (1967b), but 
are beyond the scope of this paper. Lewis (197 0) discussed estimates 
of the spectrum of counts from sectioned data, as well as pooling and 
comparison of the spectra. 



(iv) Although trend analysis was discussed in Chapter 3 of Cox 
and Lewis (1966), I don't believe it received enough emphasis. 
Also, problems such as the analysis of cyclic trends were barely 
touched upon. These and other problems in trend analysis are 
discussed in more detail later in the paper. 



(v) Finally, perhaps more emphasis could have been given to graph- 
ical methods. These were discussed in Chapter 1, but no reference was 
made to probability plotting, for example, in later chapters. See, for 
example, Wilk and Gnanadesikan (1968). 
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The field where mo*t use has been made of techniques for the 
analysis of point processes is in neurophysiological work on the signals 
occurring on nerve fibers. A summary of techniques for this type of 
analysis given in Perkel, Gerstein and Moore (1967a, 1967b) largely 
parallels Cox and Lewis (1966), with more information on the special 
problems of neurophysiology. For later work and applications, see, 
for example, Moore, Segundo, Perkel and Levitan (1970) and the 
references given in that paper. It would be impossible to try to sum- 
marize all of this work here; some of it will be touched on later, but 
there is virtually no statistical methodology given in them which is not 
given in Cox and Lewis (1966). 

Another interesting field of application is to the study of the 
occurrence of earthquakes. While strictly not a univariate point pro- 
cess, an approximation to the earthquake process as a univariate point 
process yields useful insights. For such analyses, see Vere-Jones, 
Turnovsky and Eiby (1964), Vere-Jones and Davies (1966), Vere-Jones 
(1970), and Shllen and Toksoz (1970). The discussion of Vere-Jones (1970) 
contains extensive comments on the problem of analyzing earthquake oc- 
currence data. 

Computational problems in the spectral analysis of point processes 
have been discussed by Lewis (1970). By spectral analysis here and in 
the rest of the paper we mean Bartlett's spectral analysis of the counting 
process N(t) of a point process (Bartlett, 1963; Cox and Lewis, 1966, 



10 



Chapter 5). The spectral analysis of the intervals between events in a 
point process has been greatly facilitated by the availability of the fast 
Fourier transform algorithm; see Cooley and Tukey (1965) aind Cooley, 
Lewis and Welch (197 0) for details. 

Most of the statistical methodology for the analysis of univariate 
series of events given in Cox and Lewis (1966) has been implemented in 
a computer program called SASE IV. For details, see Lewis (1966) 
and Lewis, Katcher and Weiss (1969). SASE IV is a large FORTRAN 
program with graphical output which is available from the IBM Program 
Information Department, 40 Saw Mill River Road, Hawthorne, New York 
10532, as Program No. 360 L 130001. 

In subsequent sections of this paper we discuss first two central 
problems in the statistical analysis of point processes, namely tests 
for renewal processes and tests for Poisson processes. Then techniques 
for specific non- renewal processes, such as doubly stochastic and cluster- 
ing Poisson processes are discussed. The inability to write down a 
likelihood function makes the formal analysis of non- renewal point pro- 
cesses difficult; it is only for the important doubly stochastic Poisson 
process that techniques are beginning to appear. 

The final sections deal with trend analyses, mostly of non- 
homogeneous Poisson processes; first we consider the case of monotone 
trends, then the case of cyclic trends and their relationship to the 
spectral analysis of point processes. Following this, some general 
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problems of trend analysis are considered; these include tests for particular 
rate functions and the nonhomo geneous Poisson process model per se, and 
the definition of "residuals" in the analysis of nonhomogencous Poisson 
models. 
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2. TESTS FOR RENEWAL PROCESSES 

2. 1 Markov interval processes and serial correlation coefficients, 

A natural extension of the renewal process model is to processes 
with first order dependence of intervals (Wold, 1948; Cox, 1955). The 
naturalness may be only mathematical since point processes 
of this type have not been commonly found in applications, although 
they have recently been postulated in neurophysiological contexts 
(Lampard, 1968; Walle, et al, 1969). In the process of Walloe 
et al (1969), the first order Markov interval property depends 
on the input to a neuron being a Poisson process. However, since the 
input is the superposition of an vinknown number of fairly regular neuron 
signals, the hypothesis is tenuous. A drawback in the use of the first 
order Markov interval process as an approximate model for serial 
dependence in series of events is the difficulty of obtaining analytical 
results on, for example, the spectrum of counts or the variance-time curve. 
This difficulty is closely related to the fact that there are no regenera- 
tion points in the process. 

Another problem with the model is the dearth (until recently) 
of useful bivariate distribution models. Cox (1955) used a bivariate 
exponential of the form 

f (x; X. = y) = (y) exp{-X(y)x}, (2.1) 

^i+1 
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where 



^ (y) * ^ 0^^ + X.J x) > 0, 



which had the drawback of nonlinear restrictions on the parameter. 



arises quite naturally has been discussed by Moran (1967a), Vere-Jones 
(1967), Lampard {1968), and Gaver (1970). (See also Griffiths, 1969, 
for general properties of bivariate gamma distributions. ) These bi- 
variate distributions all have positive correlations. Bivariate exponen- 
tials with negative correlation have been derived by Gaver (1972) and 
bivariate, negatively correlated intervals have been observed in a 
neurophysiological context by Walloe et al (1969). 

For statistical analysis, one of the most important theoretical 
results on bivariate exponential distributions is that of Moran (1967a), 
who showed from more general results that the serial correlation 
coefficient (or order 1) for bivariate exponential distributions is 
always greater than -0.6649; 



Another widely used model for bivariate gamma distributions which 



E{(X. - E(X))(X.^^ - E(X)} 



^ 1 - ^ ^ = -0. 6649. (2. 2) 

o 




var (X.) 
1 



Estimates of the serial correlation and tests of the renewal 



hypothesis against a first order Markov interval process are based on 
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the sample serial correlation coefficient defined as follows. 

For simplicity, assume observation starts with an event and 

ends with an event, there being n observed intervals between events 

X, , X , . . . , X with mean x » n ^ S x, . Let z. = x. - x. Then 
1 n 111 



n-1 

S z. z. , 
. , 1 1+1 
n 1*1 



1 n-1 



(2.3) 



S z. 
i*l ^ 



Before discussing theoretical results on tests for renewal pro- 



cesses based on pj^, consider its distribution for renewal processes, 
a great deal about which is now known. 

The expected value of p^^ under the renewal hypothesis, for any 
distribution of the intervals x^, is (Moran, 1967b) 



- 



1 

n-1 » 



(2. 4) 



and while its variance is known to be asymptotically n if p = 0, the 
exact variance depends on the distributions of the x. (Moran, 1967b): 



var(pj^) 



n^ - n + 1 
n^ (n-1) 



4.1 ? 

n+1 p r i ■» < n-2 

n(n-l) ^ 2 2 ^ n(n-l) * 

(S ) 



(2.5) 



Moran (1967b) obtained an approximation by replacing the expectation 
by the ratio of the expectations of the numerator and denominator. This 
result is exact for normally distributed x^'s, giving 
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var(p) 



(n-2)^ 

n^(n-l) 



( 2 . 6 ) 



For exponentially distributed x , the Moran approximation gives 

i 




n 



_7_ ^ 398 

2 3 " 4 

n n n 



(2.7) 



Moran (197 0) compared his approximation to sampling results obtained 
by Cox (1967) and a brief summary of their results on the variances 
follows. 



(a) The variance of tends to be smaller than that for the normal 
case for positive random variables with long tails. 



(b) Moran's approximation works reasonably well with gamma dis- 
tributions with shape parameter k ^ 0, n —50, or k ^ 8, n i 10. 
Outside these limits the approximation underestimates the variance 
(k = 0 is the exponential distribution). 



(c) For the exponential distribution, a partial reconstruction of 
Table 1 from Moran (197 0) using sampling results from Goodman 



and Lewis (1972) follows. 



16 



Table 2. Variance of in exponentially distributed random samples. 
Observed 

Sample variance Moran’s Difference/ 

size n (simulation) approximation Difference Observed 



20 


0. 0419 


0. 0371 


0. 0048 


. 115 


50 


0. 0182 


0. 0176 


0. 0010 


. 055 


100 


0. 00945 


0. 00935 


0. 00010 


. on 


450 


0. 00219 


0. 00221 


0. 000020 


. 009 



(d) For very skewed distributions, a better approximation to var(p^ 
is needed. This is apparent from Table 3 for random samples where 
the x^ have a Weibull distribution with shape parameter ^ ( y Wiebxall) 

A random variable with this Weibull distribution is the square of a random 
variable with a unit exponential distribution. 

Table 3. Variance of in 1/2 Weibull random samples 



Observed 



Sample 
size n 


variance 

(simulation) 


Moran* s 
approximation 


Difference 


Difference/ 

Observed 


50 


0. 01410 


0. 00058 


0. 01352 


, 96 


100 


0. 007 82 


0. 00054 


0. 00728 


. 93 


450 


0. 00201 


0. 00047 


0. 00154 


. 77 



Returning to the distribution of p^, Moran (197 0b) showed that 
n converges in distribution to a unit normal variate if the first 

four moments of the exist. Cox (1966) examined the distribution 
by synthetic sampling for the case where the x. had a normal distri- 
bution, and also distributions of the form 
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( 2 . 8 ) 



with k = 24, 8, 1, 0, - 1/2, a rectangular distribution, a double 
exponential distribution and a Cauchy distribution. The first two 
moments and coefficients of skewness and kurtosis were examined. 

For the positive random variables, the distributions were generally 
positively skewed with convergence to normality apparently fairly 
rapid. 

These results have been extended for the case of x. exponen- 
tially distributed in large scale simulations reported in Goodman and 
Lewis (1972). 

Figure 3 shows a computer plot of the estimated mean (*), 
standard deviation (x) and coefficients of skewness (+) and kurtosis 
(A) of pj as a function of n (RHO(l, N)) for n ■ 11(1) 120, 130, 140, 
150. The simulations involved at least 3, 000, 000 replications for 
each n, so that the sampling variances of the estimates are small. 
The curves have not been smoothed. Note that the skewness is positiv 
with a maximum value of about 0. 32 at n = 35, after which it starts 
back toward its asymptotic value of zero. The kurtosis is small, 
going from positive to negative at about n » 35 and then going back 
toward its asymptotic value of zero very slowly. 

The departure from normality (positive skewness) is seen 
much more clearly in Figure 4, where the computer plot gives esti- 
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mated normalized values of, from top to bottom, the quantiles at 
levels 0. 001, 0. 002, 0. 005, 0. 010, 0. 020, 0, 025, 0. 050, 0. 100, 
the estimated mean plus (n-1) \ and the quantiles at 
levels 0.900, 0.950, 0.975, 0.980, 0.990, 0.995, 0.99 8, 0.999 
These should converge to the corresponding quantiles 
of the unit normal distribution. 

Note two things: the departures from normality are relatively 

small for the inner quantiles and convergence to the normal quantiles 
is very slow. 

More detailed results from which the departure from normality 
and the slow rate of convergence can be assessed are shown in Table 4^ 
The quantiles and moments of (called RHO(l, N) in the computer) 
for n = 450 are shown in the rows marked exponential. Of these 
rows, the first row is the actual estimated value of the moment or 
quantile, the bracketed quantities in the next row are the estimated 
sampling standard deviations of the estimates (5 degrees of freedom), 

ns/ 

and the third row is the quantile minus all divided by o . This 
row illustrates the departure from normality, x^ ^ being 2. 039 
instead of 1. 96 0. 

Much more serious departures occur for the quantiles of pj^ 
from random samples with 1/2 Weibull intervals. These are given 
in the rows marked ’*1/2 Weibull'* in Table 4. Thus, not only is 
Moran’s approximation poor, as seen above, but the normal approxi- 
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mation is considerably off. The simulation results for the 1/2 Weibull 
case are discussed more fully in Goodman and Lewis (1972). 

Moran (1967b) showed that tests for independence against first order 
serially correlated intervals with exponentially distributed marginals 
based on are asymptotically most powerful. He also conjectured 
that this is true for Gamma distributed intervals. More formal results 
have been reported by Yang (197 0). 

Moran (1967b) also proposed two modifications of the serial corre- 
lation coefficient. One is obtained (see also Cox, 1955) by replacing the 
estimate of the variance in the denominator of (2. 3) by the square of the 
estimate of the mean (the mean and the standard deviation are equal for the 
exponential distribution). The other modification is due to Ogawara. 

Moran (1967b, 197 0) gives the first two moments of these test statistics 
and shows that asymptotically they involve no loss of efficiency against 
the first-order Markov interval process. However, Moran’s conjec- 
ture that the distribution of Ogawara statistic converges rapidly to the 
normal is not correct (Goodman and Lewis, 1972); in fact, its distri- 
bution is very similar to that of in the null exponential case. 

An interesting application of serial correlation statistics in 
examining neurophysiological models is given by Walloe et al (1969). 

There is very evident negative correlation between successive inter- 
vals; given the known structure of the neuronal process, first order 
Markov dependence would have been expected if the input to the 
neuron had been a Poisson process, and this was tested by checking 
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if (p.) ~ pT- Unfortunately, (personal communication) p was esti- 

mated by averaging together estimates obtained from successive sections 
of length ?0, and the bias (2. 4) accounts for a large part of the observation 
that for several experiments (p was systematically significantly larger 
than "p 2 ' jackknifed estimate (Quenouille, 1956) could have been used, 
or a correction for bias introduced. There probably still is real higher 
order dependence in this data accounted for by the fact that the input to the 
neuron is a superposition of a finite number of inputs and therefore not quite 
a Poisson process. 

2.2 Product moment score statistics 

An alternative to the serial correlation coefficient is obtained 
(Cox and Lewis, 1966, pp. 166-7) by replacing the actual interval 
values x^ by their ranks r^ or exponential scores e(r.;n) (Cox 
and Lewis, 1966, pp. 26-27), and computing a first order product 
moment statistic 

n-1 

R,(n) = S e(r.;n)X e(r ;n). (2.9) 

1 1 1+1 

A test for serial independence is then based on the null distribution of 
R^(n) under a permutation hypothesis. 

An advantage to using R^^(n) is that it controls outliers in the 
series of events; i. e. , missing points. Its distribution has been tabu- 
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lated by Lewis and Goodman (1969, 197 0) for ranks, exponential scores, 
scores from gamma populations with parameter k= -1/2 and scores 
from Weibull populations with shape parameter 1/2 (1/2 Weibull). 

The distributions for these product moment rank statistics are 
similar to those for the ordinary serial correlation coefficient from 
the equivalent parent population; positively skewed and very slow to 
converge to the asymptotic normal distribution. The distribution of the 
normalized exponential score product-moment statistic does in fact 
converge to the distribution of when the series length is n * lO, 000; 
quantiles are shown for n = 450 in Table 5 and should be compared 
with the exponential values in Table 4, 



Normalized exponential score product-moment statistic; n = 450 



Table 5. 

^0. 001 

-2.799 

’^O. 900 

1. 270 



^ 0 . 002 

-2.623 

’^O. 950 
1.656 



*0. 005 
-2. 370 

""o. 97 5 
1. 994 



"^ 0 . 010 

-2. 159 

^0. 980 

2. 096 



^ 0 . 02 0 

-1. 929 

^0. 990 

2. 395 



’^O. 025 
-1. 848 

^0. 995 
2. 669 



^0. 050 
-1. 567 

"^0. 998 

3. 009 



^0. 100 

-1. 240 

rsj 

^0. 999 
3. 245 



The sampling error in the values In Table 5 is approximately 0. 001. 

Note that the idea behind these tests is similar to that for the 
technique of random shuffling described in Perkel et al (1967a) and 
apparently commonly used in neurophysiological work. The use of the 
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product-moment statistic is more sophisticated and does not require a 
computer to do the shuffling. 



2. 3 Tests of serial independence based on the periodogram 

The first serial correlation coefficient can be used as a test 
for serial independence for alternatives other than the first order 
Markov interval process, but has certain drawbacks. In particular, 
for the more common alternatives such as clustering processes, there 
is no imperative for looking at just the first serial correlation and tests 
combining serial correlations of several orders present difficulties be- 
cause of the correlation between these statistics. 

Tests based on the periodogram were advocated by Bartlett (1954) 
and described in Cox and Lewis (1966, pp 168-17CJ and Durbin (1969). 
Denote the finite Fourier transform of the n intervals x , . . . , x by 



H ((J 
n 




1 



(Zirn) 



1/2 



n is W 

S X e ^ 
s=l « 



and the periodogram by 



U = 2irp/n, 

P 

p = 1, 2, , . , [ 1/2 n] = i- 



I (w ) = |H (w ) 
n p n p 




The test is essentially for a ”flat’^ spectrum using the distribution 
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theory for the I (w )*a. The theory states that the I (w )'s are 
n p n p 

approximately independent, exponentially distributed random variables, 
the result being exact for x^^, . . . , x^ independent and identically normally 
distributed. 

This hypothesis of a "flat" spectrum can be tested using a modi- 
fied homogeneity of variance statistic, or the cumulated periodogram 
values 



1 * 

U = 2 I <* ) / 2: I (w ) P=1,2 <-l 

' P=1 " p p.l " p 

can be tested as order statistics from a uniform distribution. Bartlett 

showed that the distribution theory is essentially independent of 

the third cumulant of the intervals between events, but it can be shown 

to be sensitive to k the fourth cumulant. 

4 

Again, results have been obtained by synthetic sampling by 
Goodman and Lewis (1972). Thedlstribution of the maximum values 
of the periodogram is shown in Table 6, the values for normally dis- 
tributed x^’s being exact, those for exponential and half Weibull being 
taken from a large simulation. The deviations from the normal case 
for exponential x^*s is small; for 1/2 Weibull x^*s, which have a 
coefficient of kurtosis of 84.720, the departure is dramatic. 

In Figure 5, for exponential varieties, we give a computer plot 
of the distribution of the modified Kolmogorov-Smirnov statistic 
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KSTWO(N) = ^^n = -s/n 



max |U -i/i| 



in the form of computer printed plots of sixteen quantiles as a f\inction 
of n, the length of the series, (n = N in the plot and the statistic is 
called KSCTWO(N)). The quantiles are a little smaller than for the 
normal case in which the upper (asymptotic) 5 percent and 1 percent 
points are 1. 35 8 and 1. 62 8, respectively. Note that for N = 140 the 
statistic is based on f = 69 U^.^'s, so that the convergence of the 
quantiles is faster than it appears to be in the figure, although not as 
fast as in the normal case. Of course, it is not known that the dis- 
tribution based on the spectrum for non-normal variates converges to 
the distribution for normal variates, but it is likely if the first four 



moments exist. 
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3. TESTS FOR POISSON PROCESSES 

We will discuss here primarily tests for the Poisson hypothesis 

against stationary alternatives, trends being discussed in Section 5. We 

assume as before that n events are observed at times t < t^ < . . . < t 

1 Z n 

in the fixed interval (Ojt^]. Intervals between events are denoted by 

x_, . . . , X , and it is convenient to denote the residual interval at the 
1 Z n 

end of the period of observation by x , = - t . 

n+1 0 n 

There are four major categories of tests. 

a) Tests based on the t^*s, conditional upon n = N(tQ), which is a 

sufficient statistic for the nuisance parameter the rate in a homo- 

geneous Pois son proces s. Since the t^*s are (conditionally) the order 
statistics from a uniform distribution, the empirical count function N(t) 
(see Figure 1) is proportional to the empirical cumulative distribution 
function for the uniformly distributed samples. The intervals x^ are 
then just the gap statistics for the sample. 

Tests based on the maximum deviation between N(t) and 
(Kolmogorov-Smirnov statistics and modifications) or other metrics 
(Anderson-Darling statistic) have well known distribution theories (see, 
for example, Durbin, 1967), but are sensitive mainly to trend departures 
from the homogeneous Pois son hypothesi s. In fact, Lewis (1965) showed 
for the special case of gamma renewal alternatives, that the test is not 
consistent. Those results can probably be extended by results on 
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empirical processes for renewal processes reported in Pyke (1972). 

The most important thing to note here is that while the distri- 
bution theory of these tests is well known because of the identity with 
the problem of testing, in a random sample of size n, a distribution 
FQ(t) against an alternative F^(t), the distribution theory under 
alternatives is completely different. In the one case F^^t^) are order 
statistics from nonuniform random samples; for point process other 
than (homogeneous or nonhomogeneous) Poisson process, both the 
uniformity and the randomness (independence) disappear. 

Note that the dispersion test for Poisson variates (Cox and 
Lewis, 1966, p. 158) is equivalent to testing the uniformity of the t^^'s 
with a chi-square goodness of fit statistic (Lewis, 1965). 

b) The gap statistics for the x^'s, i = 1, . . . , (n+1), which, under 

the Poisson hypothesis, are a random number n of independent exponen- 
tial variates, are formed from the order statistics of the x^'s as 

D . = (x .. - X ..)(n + 2 - i) (x =0, i = 1, . . . , n+1). (3. 1) 

ni (i) (i-1) 0 

These are again a random number of independent exponentials 
(Cox and Lewis, 1966, p. 26) with sum t^, and the statistics 
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t.' = 
1 



1 

S D , 
i=i 



{x,j, t + 



+ (n + 2 - )Aq (i * 1, .... n) (3. 3) 



are, conditional upon n, order statistics from a uniform (0, 1) distri- 
bution (see Cox and Lewis, 1966, p. 154 for a more complete derivation; 
the conditioning upon n and the fixed total t^ is very subtle). 

A great deal is now known about the null distributions of test 
statistics based on the gaps D ., or on the ordered x. 's, and the 
reader is referred to the excellent review articles by Pyke (1965, 1972). 
Also, against renewal alternatives, a good deal is known about asymptotic 
and, in some cases, small sample power. Alternatives are generally 
specific distributions on those in which the distribution is limited by 
specifying that the intervals have, say, distributions with monotone increas- 
ing or decreasing failure rate (Proschan and Pyke, 1967). Again see 
Pyke (1972) for a summary, and specifically Jackson (1967), Bickel 
(1969), and Bickel and Doksum (1969). 

Empirically and from a common sense point of view it is 
quite clear that against stationary alternatives^ tests for Poisson 
processes based on the L'*s are more useful than tests based on 
the L’s. This has been observed in using the SASE program. 
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in which the uniformity of the t^”s is tested using both one-sided and 



two-sided Kolmogorov-Smirnov statistics and the Cramer-von Mises 
statistic (Cox and Lewis, 1966, p. 147). , 

More information about the power of these procedures, especially 
against cluster alternatives would be valuable (Lewis, 1969} Vere-Jones, 
1970). 

c) There are also specific tests against renewal alternatives. 

The Moran test (Cox and Lewis, 1966, p. 161) is an example. 

d) A useful test can be based on the empirical spectrum of counts, 
as pointed out by Bartlett (1963). Thus, let the finite Fourier-Stieltjes 
transform of the sample fxinction N(t) be 



t 




0 



- 1/2 




dN(t) 



(3. 4) 




(3. 5) 




(3.6) 



and the periodogram 



Ij. (u) = (“^ 



0 



{A^ (w) + (w)}. 

0 ^0 
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We will refer again to the distribution theory of I (w) later. 
Noting, however, that for tot = pZir , I (to) is approximately 
exponentially distributed with mean X tt , and that for any two such 
frequencies the correlation between the periodogram points is '^l/(l+n), 
we can see that the spectral test for independence (Section 2) can be 
applied here to test for a Poisson process. 

The main drawback here is that while the spectrum of intervals 
is limited to approximately n/2 periodogram points, the count spectrum 
is not limited in such a way. Since there are roughly n/2 degrees of 
freedom available, it is a problem as to which n/2 points of the 
periodogram with frequency of the form ut^ = 27rp to use. Using more 
would invalidate, for example, use of standard distribution theory of 
Kolmogorov-Smirnov statistics to test the cumulated periodogram. 

This is a point that needs considerably more work. The tests 
are still useful in an informal way, particuarly since the shape of the 
spectrum can suggest physical reasons for departures from a Poisson 
process. In neurophysiology (Peikel et al, 1967), the tendency has 
been to use the estimated intensity function (Cox and Lewis, 1966, p. 121) 
rather than the spectrum of counts to assess departures from the Poisson 
hypothesis. This is because many of the neurophysiological processes 
causing departures are more simply expressed in time than in frequency. 
However, the distribution theory for the estimated intensity function is 
difficult. Cox (1965) has discussed the distribution theory for the Poisson 



case. 
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4. STATISTICAL ANALYSIS 

OF NON-RENEWAL POINT PROCESSES 

We consider here several non- renewal point processes which are 

of great practical importance. For the most part, however, the struc - 

tvtre is too complicated to write down a likelihood function so that 

estimation and testing is ad hoc. Examples of such analyses are cited. 

4. 1 Cluster point processes 

Cluster point processes (branching point processes) are 
important because they arise naturally in practice and have in- 
teresting mathematical properties (Lewis, 1969; Vere-Jones, 

1970). 



Briefly, a main process (usually a Poisson process) generates 
at each point a sequence of subsidiary events. In Vere-Jones (197 0) the 
main processes are initial earthquake shocks, the subsidiary or cluster 
events are aftershocks. The processes of subsidiary evens are inde- 
pendent and can, in general, be arbitrary point processes which 
terminate with probability one after a finite number of events occur. 

Two important special cases arise: 

a) The subsidiary events are generated cumulatively as a finite 

renewal case. This is known as the Bartlett-Lewis process. 
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b) The subsidiary events are generated additively, each one of the 

random number of events being independently displaced from its (main) 
generating event. The resultant subsidiary process is generally non- 
stationary and the complete process is called a Neyman-Scott cluster 
process. 



Computer failure patterns generated by this type of mechanism 
have been analyzed by Lewis (1964), and earthquakes by Vere-Jones (1970), 
Vere°Jones and Davies (1966), and Shlien and Toksbz (l970). A fairly good 
ad hoc analysis can be given for the Bartlett-Lewi s process since the 
marginal distribution of intervals is known (Lewis, 1964), as well as the 
spectrum of counts and in some cases the spectrum of intervals (Gilles 
and Lewis, 1967). The Neyman-Scott process does not yield a simple 
expression for the interval distribution and it is not yet known if the 
coefficient of variation of the intervals is greater than one, as it is 
for the Bartlett-Lewi s process. 

These cluster processes are overdispersed relative to a Poisson 
process and the variance time curve has an asymptotic form which is 
independent of the fine structure of the subsidiary processes. This is 
a help in analyzing the data; basically for large time periods the counts 
of events N(t) behave as though the subsidiary events were concentrated 
at the main, generating event, i. e. , like a bulk Poisson process. 
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Despite these points, the situation with regard to cluster processes 
is unsatisfactory; questions such as the discrimination of Bartlett-Lewis 
and Neyman-Scott processes arise in practice and are not solved (see ‘ 
the discussion in Vere-Jones (197 0)). 

When the cluster process has clusters with only one event, the pro- 
cess is equivalent to an infinite server queue. For Poisson main events 
(parameter X) the process of subsidiary events is, in equilibrium, a 
Poisson process and the delay distribution cannot be determined. 

An important special case occurs when the main process is regu- 
lar and each event is independently delayed. Appointment processes are 
very often of this type and the determination of the delay distribution from 
observation of the subsidiary events (arrivals) has been considered in 
detail by Govier and Lewis (1967) for several delay distributions. Since 
the variance-time curve has a finite limit, they have called these pro- 
cesses with controlled variability. 

4. 2 Superposed processes 

Statistical analysis of processes which are superpositions of 
point processes is important, particularly in neurophysiological contexts, 
and usually hinges on questions of estimating the number of contributing 
processes or, when this is known, identifying the structure of the com- 



ponent processes. 
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These are difficult questions whose solution has not progressed 
much beyond the basic work of Cox and Smith reported in Cox and Lewis 
(1966, Chapter 8). Without specific assumptions; i. e. , that the com- 
ponent processes are renewal processes with, say, gamma distributed 
intervals, very little can be done to estimate the number of processes. 

The problem of Walloe et al (1969) is of this kind, the nature of the 
input processes and the neural mechanism being well known, the main 
question being how many inputs impinged on the neuron. 

Identifying the component processes is again, difficult, without 
specific assumptions. Work such as that of Ross (197 0) on identifying 
the interval distribution in a known number k of superposed processes 
is rather technical; note that the variance time curve and spectra of 
counts are additive; e, g. , the spectrum of counts of the superposed 
process is k times the spectrum of the individual processes. Thus, since 
the spectrum of counts of a renewal process is simply related to the Laplace 
transform of the probability density of the intervals, f*(s), 

^ f«(iw) f*(-iw) 

X X 

where |i is the mean of x, the spectrum of the superposed process 
g^((j) = kg^((j) and it is natural, and probably efficient, to use the 
estimated spectrum or the estimated covariance density to estimate 
f (x). Convergence follows from results of Brillinger (1972). 
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Nonhomogeneous superpositions also occur (Blumenthal, Green- 
wood, and Herbach (1971)), but as in the homogeneous case, when the 
number of contributing processes is large, the resultant process over 

I 

realistic periods of observation is almost indistinguishable from a 
(homogeneous or nonhomogeneous) Poisson process. 

4. 3 Doubly stochastic point processes 

Let A{t) be a real valued, nondecreasing stochastic process. 

A doubly stochastic point process is a generalization of the nonhomo- 
geneous Poisson process in which the integrated intensity function is 
replaced by A{t). Thus, given a realization of A(t), the point process 
is a nonhomogeneous Poisson process. Generally Alt) is differentiable 
and X (t) = A^[t) might be a stationary stochastic process, a determin- 
istic trend or a combination of both. The process was introduced by 
Cox (1955); see Bartlett (1966, p. 325), Cox and Lewis (1966, p. 179), 
and Grandell (197 0) for details. Gaver (1963) considered the case where 
X (t) changes level at random times and called the process a random 
hazard process. 

The doubly stochastic mechanism in a point process is very 
realistic and probably quite common. Thus computer failure processes 
depend to some extent on temperature, humidity, etc. Unfortunately, 
analytic properties of the process when X (t) has a stochastic element 
are very difficult to derive; see Lawrance (197 2). If X (t) is a stationary 
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stochastic process with mean variance o and autocorrelation 

K 

function 



OO -iMX 

p(T) = l/?7r / p(r) e dF(u), 

- OO 



( 4 - 2 ) 



where F(cj) is the integrated spectrum and f(y) =F*(cj) when it exists, 
then several useful results can be obtained (Cox and Lewis, 1966, 
p. 179-183). In particular, the covariance density \^(t) of the doubly 
stochastic Poisson process is 



V+(t) = p(t). 



(4. 3) 



and the spectrum of counts 



g^(u) = X/tt + 2f(U)). 



(4. 4) 



Also, in general, the index of dispersion is 



2 ^ 

I(oo) = 1+2 > / p(u) du, 

0 



(4. 5) 



SO that the process is overdispersed relative to the Poisson process. 

There has been a great deal of interest in the doubly stochastic 
Poisson process in physics, optics, and engineering, and some work 
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on the estimation problems. 

a) In physics photon emission is a well understood physical process, 

and it is known that the process modulating the Poisson emission is, for 
instance, an Ornstein-Uhlenbeck process. The parameters of this 
process are physical parameters which it is of interest to estimate, 
and attempts have been made (e. g. , Pusey, 1971, Jakeman, Pyke and 
Swain, 1971, and Koppel, 1971) to do this via the spectrum and the 
relationship (4. 4), Again, laser light is deliberately focused on a physical 
system and the resultant intensity fluctuations in the scattered light re- 
flects rates of molecular motions and interactions in the system. 

In most cases the number of photon counts in these experiments 
is very large and it is convenient, and very often necessary to cumulate 
counts. There are problems of determining the best sampling interval, 
a problem which is often made simpler by detailed knowledge of the 
modulating process. 

Perhaps because of the large amount of data involved and the 
fact that the X (t) process is fairly well understood, physicists rarely 
bother about details such as the efficiency or optimality of estimation 
procedures. Moreover, the spectral estimation is very often automated 



and done digitally. 
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b) In engineering, two areas are involved. One is biomedical 

engineering in which, for example, radioactive substances are injected 
into the blood stream and the counts of the radioactive emission are used 

I 

to estimate a decay function which is related to a physical process of 
interest (Snyder, 19*71). 

A second area is optics, where modulated light is used to trans- 
mit information (Reiffen and Scherman, 1963; Bar-David^ 1969; Karp 
and Clark^ 1970; Clark and Hoversten^ 1970), Here \ (t) is very often 
a signal changing levels, a possible problem being the discrimination 
of several X .(t)*s from photon counts of the noise. Very often a 
stationary **noise'* element is also present; see Bedard (1966) for the 
physical considerations. 

In all of this engineering literature, there is much concern 
with optimality, in some sense, and procedures are based on likelihood 
ratios, Bayesian posterior statistics and maximum likelihood detectors. 

It is a difficult literature to penetrate and my overall impression is that 
there are many hidden assumptions involved, one being that samples 
are large, the other a normality assumption which may be quite incorrect. 
Moreover, most explicit results are for very simple situations such as 
mixed Poisson processes of one type or another, or rates (Reiffen and Sherman, 
1963; Snyder, 1972) changing at known times. The latter problem is very simple 
and straightforward, e specially when compared to the case of unknown change 
points which gives a true doubly stochastic Poisson process and an 
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inference problem similar to the difficult change point problem. 

c) Grandell (1971, 1972) has considered inference in doubly stochastic 

processes from the viewpoint of mathematical statistics, his motivation 
apparently being primarily problems in actuarial work. His approach is 
quite different from that in the engineering literature, cited above, and 
he has considered, for example, optimum estimates of E{X.(t)} based 
on linear combinations of the data. The problems are related to the 
general problem of curve estimation from random data and, while 
Grandell (197 2) gives specific examples, much remains to be done. 
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5. TREND ANALYSIS 

IN NON-HOMOGENEOUS POISSON PROCESSES 

We discuss now the analysis of trends in non-homogeneous Poisson 

processes, extending the results of Cox and Lewis (1966, Ch. 3). Inltlallyi 

we discuss results based on specific parametric models for the rate function 

X(t) of the non-homogeneous Poisson process. These results are based on 

the fact that the likelihood fvinctlon for n observa.tlons In the fixed period 

(0, t ] at times t,<t_< <t Is 

*0 12 n 

n ^0 

L(tj, . . . , tn* n; p) = 1^1 \(ti; e ) exp { - X(u;e) du}, (5.1) 

where ^ denotes the vector of parameters In the model. Moreover, given 
that n events occur in (0,t^], the times to events t^ are the order statis- 
tics from a random sample from the probability density function 



f(t;6) 



X (t;^ 

ox (u;0)du 



X (t;Q) 
A(tQ;S) 



0 < t < t 



0 ’ 



(5.2) 



and the conditional likelihood Is 



L(tj, . . . , t^;n;0) = n| X(t.;9) / (A(tQ;0)}’ 



(5.3) 



Later In the section we discuss procedures for examining the adequacy 
of t)^e model for X(t) and adequacy of the non-homogene ous Poisson process 
model Itself. 

5. 1 Monotone and evolutionary trends 

The estimate of X(t) when there Is no trend present 1. e. X(t) = X, Is 

n/to- 



45 



In Cox and Lewis (1 966,. Chapte r ^ the exponential linear trend 



X(t) = exp t} = X exp {fi t} 



(5.4) 



was discussed and It was shown^ using (5. 1), that (n,2t^) were a sat of 



sufficiant statistics for (or, ^). A test for =0 against ^0, whera a is 



is an optimum (conditional) test. The maximum likelihood estimate of 
was also given implicitly, but no small sample properties of the estimate 
are known. 

The test for /5 = 0 was applied in Section 1 to analyze a sequence of ar- 
rivals at an Intensive care unit (Table 1) and a sectioned analysis indicated 
that the trend was not monotone, although same overall Increase In the rate 
was present. 



The maxlmixm likelihood estimates of X, a,/J are given as the solution of the 
equations 



a nuisance parameter, is based on the distribution of lit./n, given n. This 



An exponential quadratic trend uses the model 




(5.5) 



for which we get 




(5.6) 




(5.7) 




n 



(5.8) 
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^.2 . 0 2 - - 2, , 

Stj Jq u ®xp ifin + yn ) du 

-r “ 1 = (5.9) 

/Qexp(^u+^u ) du 

2 ' 

and it is clear that (n, I^t., I^t^ ) are a set of sufficient statistics for 

K P* V). 

There are several interesting open questions here. 

What are the small sample properties of the estimators and what are 
their large sample properties? Note that the usual theory for maximum 
likelihood estimates is not directly applicable because of the random sample 
size. 

What effect does the quadratic term (1. e. \ / 0) have on the estimates 
of p in the model (5.4)? 

For the arrival data the following estimates were obtained. 

Linear model: X = 0. 6342{ j? * 0001 427j ^ , 

Quadratic model: X « 0.4792} p = 0. 0009788: ‘Y=-0. 4481x10 

Note the large difference in P in the two cases. 

Another problem of interest is to test for the quadratic term in the 

trend. It is clear from Table 1 that either a quadratic team or a long cycle 

is necessary to adequately model the annual data from Section 1. 

A formal tast csui be based on the idea that» for any given n and 2^t^ 

are a set of sufficient statistics for or,/3. Therefore a test for \ is based 

2 

on the statistic Zt^ and Its conditional distribution, given n and Ft^. This 

2 

distribution is difficult to obtain analytically, Z)tj being the square of the 
distance to the sample point, which is constrained to lie in the (n-1 )-dimen- 
slonal hyperplane defined by Zt^ = C. 
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For large samples this distribution can be obtained from the fact 
2 

that 2tj/n and 2t^/n, conditioned on n, are jointly normally distributed 
for large n, and the following exact moment results; 



= E{ Ftj/n} = t^/Z; o-^ = var { S tj/n} = t^ (1 2n) S 

H -2 = E{Zt^/n}= t^/3; cr^ = var { Et^/n) * 4t^ (45n) 

/ T n/15 - qAq 

p= corr{— 1, — * } = — r = 0. 9o8. 

n n 4 

2 

Thus, using normal theory results we test with Et^/n having a normal 
distribution with mean and deviation 



“■l / 2 \ 

|x = jx,+pji (_1 - 1^,.) 






2 n 



= >0'^ +‘o'— ’ - ^ 

n 2 



<r =<r^ (1-p )*. 



(5. 10) 
(5.11) 



For the complete set of arrival data Et^/n = 954. 25, giving p. = 1, 1 87, 783, 
and O’ = 6, 530, Z tf/n = 1, 161, 565, giving (Zt^/n-p)/o- =-4.015 and this 
is clearly a highly significant rejection of the hypothesis that y = 0. 

Boswell (19 66) and Boswell and Bruiik (1969) have considered tests 
of the hypothesis that X(t) is not constant but is non-decreasing. Using a 
likelihood ratio criterion, conditional on the fixed value of N(tQ)»n, he 
found the test statistic 

where 
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iiftl I?K^n - V’ 

the null hypothesis being rejected for large values of the statistic. 

The statistic (5. 12) is not simple to compute, but Boswell gave an ' 
iterative procedure for the computation, and same results on the limiting 
distribution of the statistic. 

It would be of interest to compare the power of this test against the 
power of the test based on Zt^/n for the exponential linear trend. Some 
results have been obtained by Mr. Ian White of the University of Edin- 
burgh. 

5. 2 Cyclic trends of fixed frequency 

Cyclic time trends (as opposed to cycles on serial number) occur 
frequently in point processes but were treated only as an exercise in Cox 8c 
Lewis ( 1966 ). An example of such a series is given in Forrest (1950), who 
was investigating thunderstorm severity in Great Britain and its effect on 
power lines. He found a tendency for thunderstorms to occur in the morning 
(time of day effect) and of course a very strong time of year (seasonal)effect. 

The arrival data discussed in Section 1 might be expected to have 
such effects, although a time of day effect is by no means evident since there 
is only about one arrival every day and a half. In Figure 6 we show a 
collapsed plot of the nvimbers of arrivals in the successive hourly periods 
of all days of observation. The plot should be compared to Figure 1 . 8a in 
Cox 8c Lewis (1966) for the arrivals from 4 February 1963 to 1 8 March 1963. 
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Figure 6. Arrival of patients at an intensive care unit. 

Collapsed hourly plot of number of arrivals to investigate ”time-of-day " effect. 
Second period of observation: 19 March 1964 to 8 February 1968. Total number 

of arrivals (including ties): 1216. Observation time t^ *= 1420 days; m = 0. 856. 




DAY OF WEEK 



Figure 7. Arrival of patients at an intensive care unit. 

Collapsed daily plot or arrivals. All arrivals including ties. Solid line 
is 2nd period, 19 March 1964 to 8 February 1968 (1216 arrivals in 1420 
days). Dashed line is complete record, 4 February 1963 to 8 February 
1968 (1467 in 1829 days). . Z ^ (n. - n)^/(n) = 12. 93, Note that 
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The plot for the complete period (4February 1 963 - 6 February 1968) in 
Figure 6 is similar in shape to the plot for the earlier period and no formal 
statistical test is needed to conclude that there is a real time-of-day effect;. 

As a model for the rate \(t) in a non*»homogeneous Poisson process 
with fixed frequency u)Q=2tr/TQ we take 

\(t) = exp {or + k sin w t+k cos u t} (5.14) 

S U C (J 

= X. exp {k sin(UQt + 6)}, (5.15) 

V^) 

2 2^ -1 ^c 

where k=^^+k^}^, 6 = tan ( — ) , X. = exp(or)x I^(k), 



and Iq(^) Is a modified Bessel function of the first kind of zero order. This 

reparameterization is used because 
T 

/ ° exp {k sin ( -^*^ + e ) } dt = IqOc) , (5.16) 

0 



so that if tp, the total period of observation is a multiple p of T^ (here one 
day), we have 






\(u)du = X t^ = XpT^. 



(5.17) 



In effect we have separated out multiplicative ly a linear growth from a 
cyclic effect which takes place only within periods of length T^. 

The reason for using the form 5. 14 rather than, say, 

\(t) = ar + k sin (u^t +0 ) (5.18) 

is that X(t) must be positive, and to achieve this with (5. 18) requires k<ar. 
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For k « Of the two models are essentially the same. In 
addition, the rate (5. 14) gives simple results based on sufficient 
statistics, this being closely connected with the fact that the density 
(5. 2) belongs to the exponential family of density functions. 

An additional reason for using (5, 14) is that such nonsymmetrical 
cycles are probably more natural with series of events, possibly because 
of the positivity of the rate function* Professor J* W. Tukey has suggested 
a model which is (5* 18) squared; this could be useful with data in which there 
are a large number of arrivals in each period, so that regression techniques 
using the square root of the numbers of arrivals in fixed intervals can be 
used (Cox and Lewis, 1966, Chapt. 3), This model has also been used by 
Fisher (1953). 

In Lewis (197 0) the following results were given for the non-homo- 
geneous Poisson process with the rate \ (t) given in (5. 14). 

Using (5. 3) we find that the observations enter into the likelihood 
function only through 



n, 



t ' 0 
0 



n 

S cos co^t., 
^1 0 



‘o " 



n 

S sin CO t. , 
i=l 0 ^ 



(5. 19) 



and these are a s«t of sufficient statistics for the parameters a, 

and k in (5. 14). Note that (co.) and B (to ) are the components 

of the periodogram (3. 6). 

Maximum likelihood estimates of X , 9, k are 
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9 = t an 




(cOq) 



(cOq) 



(5. 20) 



X. 




(5. 21) 



and k is the solution of the equation 



(l/n){A^ (CO ) + (CO )} 
0 0 



I (k) 

= »lj(k), 

io(k) 



(5.22) 



where Ij^(k) is the derivative of and i(;(k) increases monotonically 

from 0 to 1 as k goes from 0 to oo. This later fact allows one to 
use the Neyman-Pear son lemma to show that, conditional on the observed 
value n of N(t^), the most powerful test of k = 0 against k 0 is 
based on the statistic 



(•o/n) 



t,2 , , 

B (co.) 

‘o “ 



)= (iTt /n) I ). 

0 



(5.23) 



Since, for k = 0, 



21 (co„) 

‘o " 



has asymptotically 



(Cox and Lewis, 



1966, Chapter 5) an exponential distribution, the test for k = 0 reduces 
to accepting the hypothesis at, say, a 5% level if 



21^ (cOq) = (6/iT)X(tQ/n)= {6 /tt)X m . 



(5-24) 



The factor 6 arises because prob \\ = 6} = . 95, where is 
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random variable with a chi-square distribution of two degrees of freedom. 

This test applied in Lewis (197 0) to the first section of the arrival 
data gave a strong indication of the presence of the cycle at a period of 
one day. 

For the complete record of arrivals at an intensive care unit dis- 
cussed in Section 1 we get for the periodogram at p = 1829 or = 2 tt, 
which is the day frequency, since t^ = 1829 and 0 )^ = 

ZI (co ) = 27, 094 and 6 m/ 7 T = 1. 52. 

''0 

This is, not surprisingly, highly significant. For the first 409 
arrivals the results are (Lewis, 1970) 

21^ =: 5. 0 and 6m/7r = 1, 10, 

Thus the periodogram component is increasing roughly in proportion to 
n, as it should for a true cyclic component. 

When Lewis (197 0) was written, the connection of this model for 
a cycle of fixed frequency with tests for directionality on the circle was 
not realized. In fact, the conditional test (5. 24) is equivalent to testing 
that n observations on a circle (here a 24-hour clock) have the von Mises 
or circular normal distribution (Gumbel, Greenwood and Durand, 1953) 
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f(i) = 



exp fk sin(i + 9) } 
2 IT I^(k) 



(0 ^ 2tt) . (5. 25) 



For k = 0, this is a uniform distribution; otherwise, it has a mode at 
9, the vector in the direction 9 being called the modal vector. Green- 
wood and Durand (1955) obtained the distribution of the square of the 
resultant 



Z ^ 2 

R = ( S cos i.) + 

i=l 1 



2 

( Z sin i .) 
i=l 1 



(5. 26) 



of n such vectors when k = 0, generalizing earlier work of Pearson 
on the problem of random walks on the circle. The formal analogy of 
(5. 26) with the periodogram (3. 5) should be clear. Watson and 
Williams (1956) generalized these distributional results, using results 
for sufficient statistics, and found the conditional p. d. f. of the quantity 
on the left of equation (5. 22), which we denote by r, as 



f(n^^^r) = r I^(kr) / (j^(u) jQ(ru)du/ (lQ(k) (5.27) 



oo 

0 



where Jq( * ) ordinary Bessel function of zero order. 

It is not all apparent that (5. 27) is more useful for computation 
than the generating function, which is given, for example, for k = 0 
in Cox and Lrewis (1966, Chapter 5) in the discussion of the distribution 
of the periodogram at one point where is of the form 



55 



= 2-n-p, for p integer. However, Stephens (1969a, 1969b) has 
used (5. 27) to tabulate, among other things, the power of the test that 
k = 0 against alternatives k > 0. The most complete discussion of these tests 
for the von-Mises distribution is given in Watson and Williams (l956). The 
function 4^(k) is tabulated by Gumbel, Greenwood, and Durand (l953). Stephens 
(1969b) has also discussed tests for 0 = 9^, and joint tests for k and 9. 

It is clear that problems involving cycles at two or three fixed 
frequences, e. g. , to , to , to , will arise in analyzing series of events. 

For example, in the problem of analyzing the arrivals at an intensive 
care unit, a time of week effect and a time of year effect are distinct 
possibilities. Surprisingly enough, there does not seem to be a strong 
time of week effect in the data, but there is not space here to go into 
this. 

Formal tests for more than one cycle follow from results in 
Cox and Lewis (1966, Chapter 5) to the effect that at two different 
frequencies and to^, both of which when multiplied by t^ are Ztt 

times an integer, the correlation between the periodogram values 
(cJq) and I (a)^) is approximately (1 + n) \ Thus, a test for a 
time of week effect based on the conditional distribution of 

given the values of N(t ) = n and I reduces in effect to an inde- 

0 tp 0 

pendent test at based on (5. 23). 

For the arrival data, considered over a period of 1456 days, 
we get for p = 208 (to * 27 t/ 7), or a period of one week, the value 
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(to ) (to ) 

2 P ^0 P 

21 (u ) = - {— + — } = 0.467. 

‘o P " ‘o ‘o 

Thus there is no formal indication of a time«of«week effect. 



5, 3 The spectrum of counts and cycles of unknown frequency 

In the previous subsection, we showed the connection between 
tests for a cyclic component at a known frequency in a nonhomogeneous 
Poisson process and the periodogram, the periodogram being the basis 
for estimation of the (second order) spectrum of counts g^(co). 
Clearly, one might want to look for cycles at unknown frequency and 
this will, intuitively, be based on the spectrum of counts. More pre- 
cisely, it will be based on the periodogram 



(w) = (1/irtp) (A^ (O)) + (to)}. (5.28) 

0 0 0 

The analogous problem in ordinary time -series analysis is the 
classical problem of hidden periodicities, discussed at length by Hannan 
(197 0, p. 463). This problem has not been tackled in point processes, 
and is more difficult for two reasons. First, the periodogram points 
are not quite uncorrelated, as we saw in the previous section, and also 
the spectrum is, in theory, not limited in the frequency domain. (There 
will always be some band limiting due to jitter; see Lewis, 197 0. ) 

Thus, it is a problem how to use the distribution theory for the spectrum 
given in Bartlett (1963) and Cox and Lewis (1966, Chapter 5) and how to 
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pick relevant parts of the spectrum from which to estimate the unknown 
frequency. 

As an example, we give in Figure 8 the spectrum of the arrival 

data of Section 1. The peak for the day effect at p = 1829 (co = Ztt) is 

P 

evident, but it would be difficult to read into this spectrum anything but 
a conclusion that it is a nonhomogeneous Poisson process with a single 
fixed cycle (day-effect). There is possibly a harmonic at p = 3658 
and some inhibition at low frequencies, but all the points (except for 
p = 1829) are within 1 percent confidence bands for individual values, 
and would be within bands for the maximum of the periodogram values 
(Bartlett, 1963). 

The spectrum of counts is generally a useful tool, and even more 
particularly so for non-Poissonprocesses , but has found very limited 
acceptance amongst applied workers. This is partly due to confusion 
with the spectr\im of intervals, with which neurophysiologists, for example, 
are much more familiar. 

Lewis (197 0) has given a heuristic justification for the spectrum 
of counts and discussed computation and smoothing. Bartlett (1967) has 
discussed finer points in this type of analysis. 

Finally, Brillinger (1972) has given a general theory for spectral 
analysis of point processes, and has defined higher order spectra 
(cumulant spectra of order k) for point processes. These are gener- 
alizations of the cumulant spectra of order k of continuous time series. 
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p(oip = 27rp/to) OR PERIOD T IN DAYS 



Figure 8, Arrival of patients at an intensive care unit. 

Spectrum of counts for the complete record: 4 February 1963 to 8 Febraury 1968. 

n « 1458 arrivals in a period of t^ = 1829 days. 60 point uniform smoothing. 
Bands are 0. 95 and 0, 99 confidence regions for individual values, m * 0. 797 2. 




SERIAL NUMBER OF MIDDLE ARRIVAL IN GROUP OF 20 INTER-ARRIVAL TIMES 



Figure 9. Arrival of patients at an intensive care unit. 
Averages of successive groups of 20 inter-arrival times after detrending with 
rt ^ 

T = \ exp { a + ^ + k sin ( 2 it u + 0 } du, where a, /3 , k, 0 are maximum likeli- 
hood ^estimates, n *= 1458 arrivals in 1829 days from 4 February 1963 to 
8 February 1968, Overall average after detrending: 1. 0088, Homogeneity 

of variance statistic = 309. 07 (x ^ mean 71, o-*= 11.92), 
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Brillinger's work is based on a general spectral representation for N(t), 
the counting process of the point process, and provides an asymptotic 
theory for the spectral estimates of order 2 (periodogram) considered 
in this paper. Brillinger’s paper is far too extensive to do justice to 
here. The problem of whether cumulant spectra of order higher than 2 
will be useful in practice remains open and should be explored. 
Brillinger’s spectral decomposition may provide an answer to whether 
spectral analysis can be useful as a representation in real systems; i. e., 
neurons. 



5. 4 Mixed models 

The arrival data considered in Section 1 has been shown to have 
a long term trend which can be represented by the exponential quadratic 
function (5. 5) and a strong day cycle. A combined model for this data 
could be a nonhomogeneous Poisson process with rate X (t) of the form 

X (t) = exp (or + pt + yt + k sin(o)Qt + 9)}. (5.29) 

t 

This is difficult to handle because A(t) = / ^ (u)du is not tract- 

0 

able. However, in the context of the arrival data the evolutionary trend 
changes little within a cycle and an approximation in which the linear 
and quadratic terms are assumed constant over the period of each cycle 
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yields a tractable model. The estimates of k and 0 are then simply (5, 20) 
and (5. 22), y and p are estimated from (5. 8) and (5. 9) and X is 
given by (5. 7) with lQ(k) multiplying the term in the denominator. 

For the arrival data, we obtain the following values: 

X = 0.47 92; 

k = 0. 654; 

9 = 3.519; 
p = 0. 0009788; 
y = -0. 4481 X lO'^. 

Clearly, mixed models such as (5. 29) will be needed in many 
cases for the analysis of real data. 



5. 5 Residual analysis of point processes 

Many observed point processes have rate functions which are 
clearly not easily representable by simple, tractable models such as 
those discussed in the previous section. For example, it is common 
in observing arrivals of jobs in a computing center to find (roughly) 
the rate increasing sharply throughout the morning, dipping around 
Ixinch time, picking up slightly in the afternoon, and then dipping off 
sharply at night. The situation is fairly simple in this example, since 
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many days observations are available, showing similar variations, 
and careful pooling gives a smoothed tabulated estimate of the rate. 

This rate is needed to model scheduling algorithms for the computer. 

In general, the problem of smoothing the point process to obtain 
an estimate of X. (t) is, in my opinion, open, despite the work on doubly 
stochastic Poisson processes cited in Section 4, which should be applic- 
able. It is generally connected with the statistical problem of curve 
smoothing, and will not be considered further. 

Both in the case where \ (t) is obtained by smoothing or from a 
specific model, say exp {a + pt}, it will often be necessary to test the 
adequacy of both the trend model and/or the assumption of an underlying 
nonhomogen'eous Poisson process. 

For example, an intensive care unit is a service facility, and to 
provide adequate service, one must know something about the arrival 
process; it might be very regular if all arrivals are from surgical 
operations, or exhibit clustering if all arrivals are from the emergency 
room. The form of the underlying process will also affect the test for 
cyclic trend given in (5. 23), the attempt to estimate the spectrum and 
filter out the spike usually being a selfdef eating process (Bartlett, 1967; 
Lewis, 197 0). 

Both formal and informal methods, usually graphical, are re- 
quired for this question, analagous to the similar problem in regression 
analysis which is usually approached by examining the residuals after 
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estimating parameters for the mean value in the model. 

While keeping in mind that there are many ways to approach 
this in different situations, the following proposal seems to offer a 
systematic approach to examining such questions (Lewis, 197 0). It is 
well known that on the transformed time scale 

t 

T = / X (u)du, (5. 30) 

0 

a nonhomogeneous Poisson process with rate function X (t) becomes a 
homogeneous Poisson process of rate 1. Thus, we proposed to examine 
the residual process where 

t- ^ 

T. » r ^ X (u)du, (5. 31) 

1 0 

and X (t) is some estimate of X (t). 

There are a number of ways of looking at this transform. 



1) If X (t) is very smooth over periods compared to the mean time 

between intervals, we can write 



X.= T. 



- T 



i-1 






(5.32) 



Then a logarithmic transformation reduces this to standard regression 
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models, and the exponential trend functions discussed above are appealing 
since we have (approximately) standard linear regression methods. 

These are considered in Cox and Lewis (1966, Chapter 3). 

Of course, in data such as the arrival data of Section 1, this 
approximation is not possible, and the integral transformation must 
be performed, perhaps nvimerically. 

The real problem comes, of course, when X (t) cannot be 
assumed to be smooth over suitably short intervals and no simple 
parametric form for X (t) is available. 

2) If the analysis is based on the conditional distribution of the 

t^'s, given N(t^) = n, then the transformation (5. 31) is seen from (5. 2) 
to be just proportional to the probability integral transform for the density 
X (t)/A(t) with estimated parameters. 

3) There is also a possibility of looking at (5. 31) as a filter in 
the spectral domain. 

Formal methods for analyzing the are difficult and will be 

discussed elsewhere. This is simplest for X (t;9) when sufficient 
statistics for 9 are available. The probability structure of the 
process, given the sufficient statistics, is independent of Of 

course, in the simplest case of X (t) = X , with X * well 



64 



known that the X V s are the gap statistics for independent uniform 
random variables and are asymptotically exponentially distributed 
and independent (Pyke, 19?2)« For small samples, they will gener- 
ally be correlated random variables. 

Informal, graphical methods are useful and are discussed in 
the context of the arrival data given in Section 1. The first section 
of this data was transformed, using (5. 31) and an exponential linear 
trend, and discussed in Lewis (197 0). No test for Poisson processes 
performed in the SASE IV program gave significant indication of 
departure. The spectrum of the transformed data is shown in Lewis 
(1970). 

Clearly the (exponential) linear trend plus day cycle is not adequate 
for the complete set of data, as we saw in Section 1, We examine its adequacy 
further here. Figure 9 shows a plot of the average of successive groups of 
20 after detrending with the (exponential) linear trend plus day cycle. The 

model for X(t) is clearly inadequate. Figure 10 shows the spectrum of 
counts after the linear transformation. There is no spike corresponding 
to the one day period and no Indication of departure from a Poisson process. 
Similarly, the tests discussed in Section 3 for Poisson processes do not give 
any great indication of departures. However, the estimated asymptotic slope 
of the variance-time curve Is approximately 8, and since this quantity is re- 
lated to the rate of change of the low frequency end of the spectrum, there is 
again an indication of inadequacy of the form of X(t) or the hypothesis of a 



nonhomogeneous Poisson process. 
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A quadratic term was added to the rate X. (t), giving the trans- 
formation 

^i »A A 2 •^.1 

T . = J exp for + pu + 'y'U + k sin(2Tr u + 9) }du, (5. 33) 

^ 0 

where the values of the estimated parameters are given in Section 5* 4. 
The plot corresponding to Figure 9 is shown for the (exponential) 
quadratic detrending in Figure 11. This and other results give a much 
clearer picture of a possible cycle or quasi-cycle in the data with period 
of about a year and a half. The example is discussed in more detail in 
Lewis (1971), where local smoothing is used to estimate X (t) for this 
quasi-cycle and test the nonhomogeneous Poisson hypothesis. 
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p(cup=2Trp/fQ) OR PERIOD T IN DAYS 

Figure 10, Arrival of patients at an intensive care unit. 
Spectrum of counts for complete record after detrending with 

AAA A ^ ^ A A 

T = \ exp /5u + ksin{27ru + 0)} du, where a, /3 , k, 0 are maximum 

likelihood estimates, n = 1458 arrivals in 18?9 days from 4 February 1963 to 
8 February 1968. Bands are approximately 0.95 and 0.99 confidence regions 
for individual values. 




SERIAL NUMBER OF MIDDLE ARRIVAL IN GROUP OF 20 INTER -ARRIVAL TIMES 



Figure 11. Arrival of patients at an intensive care unit. 
Averages of successive groups of 20 inter-arrival times after detrending with 

7A A y.^AAAA 



k, 



are maxi- 



T = \ exp {a + ^ u + y u^ + k sin { 27 t u + 0) } du, where o', /? , y , 

mumS.ikelihood estimates, n = 1458 arrivals in t^ = 1829 days from 4 February 
1963 to 8 February 1968. Overall average after detrending: 

Homogeneity of variance statistic: ( mean 71, cr » 11.92). 
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6. CONCLUSIONS 

Much further work is needed in the statistical analysis of series 
of events. Particular areas that come to mind are the analogues of 
the procedures in Section 5 for modulated renewal processes. Likeli- 
hood functions can be written down, but results are fairly difficult to 
obtain from the expressions (see Cox (197 2) for details). Both for these 
processes and the nonhomogeneous Poisson process there is a great 
deal of work to be done to justified estimation and testing procedures based 
on random sample sizes. Some of these justifications have been given by 
Brown ( 197 2). 

Finally, it is hoped that new methods will be developed which will 
allow some of the nonrenewal point processes to be analyzed in a more 
systematic way, 
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